# Appendix A Balance Testing

# MISSING the original code for balance testing figures. Here is some default code that will show no balance
# issues exist (I paid for the data, so it is of the highest quality).
# 1 variable (gender in study 2b) comes up as ever so slightly unbalanced (p=0.035), but having 1 in 48 variables come up
# as unbalanced in consistent with statistical noise. A quick eye-check at the figures in the main text will
# show that whatever imbalance does exist must be extremely slight.

library(cobalt)

txdata <- read.csv("./tex.csv")
library(dplyr)
tx_pos <- txdata %>% filter(treatment < 2)
tx_neg <- txdata %>% filter(treatment > 1) %>% mutate(treatment = treatment - 2)

# p-values showing balance for each study's variables 1-by-1:

# study 1a
vars <- c("Age","Gender","Education","Income","Born.in.U.S.","Partisanship",
          "Travel","SDO","txknw","Family","Duration..in.seconds.","Mexican")

res <- data.frame(var = vars, p = NA_real_)

for(i in seq_along(vars)) {
  v <- vars[i]
  x <- tx_pos[[v]]
  
  # numeric vs categorical branch
  if(is.numeric(x)) {
    res$p[i] <- t.test(x ~ tx_pos$treatment)$p.value
  } else {
    tbl <- table(x, tx_pos$treatment)
    res$p[i] <- chisq.test(tbl)$p.value
  }
}

print(res)

sum(res$p < .05, na.rm = TRUE)

# study 2a
vars <- c("Age","Gender","Education","Income","Born.in.U.S.","Partisanship",
          "Travel","SDO","txknw","Family","Duration..in.seconds.","Mexican")

res <- data.frame(var = vars, p = NA_real_)

for(i in seq_along(vars)) {
  v <- vars[i]
  x <- tx_neg[[v]]
  
  # numeric vs categorical branch
  if(is.numeric(x)) {
    res$p[i] <- t.test(x ~ tx_neg$treatment)$p.value
  } else {
    tbl <- table(x, tx_neg$treatment)
    res$p[i] <- chisq.test(tbl)$p.value
  }
}

print(res)

sum(res$p < .05, na.rm = TRUE)

# study 1b

mindata <- read.csv("./min.csv")
library(dplyr)
library(stargazer)
min_pos <- mindata %>% filter(treat < 3) %>% mutate(treat = treat - 1)
min_neg <- mindata %>% filter(treat > 2) %>% mutate(treat = treat - 3)

vars <- c("Age","Gender","education","income","us_born","partisanship",
          "Travel","SDO","minkn","Family","Duration..in.seconds.","white")

res <- data.frame(var = vars, p = NA_real_)

for(i in seq_along(vars)) {
  v <- vars[i]
  x <- min_pos[[v]]
  
  # numeric vs categorical branch
  if(is.numeric(x)) {
    res$p[i] <- t.test(x ~ min_pos$treat)$p.value
  } else {
    tbl <- table(x, min_pos$treat)
    res$p[i] <- chisq.test(tbl)$p.value
  }
}

print(res)

# study 2b

vars <- c("Age","Gender","education","income","us_born","partisanship",
          "Travel","SDO","minkn","Family","Duration..in.seconds.","white")

res <- data.frame(var = vars, p = NA_real_)

for(i in seq_along(vars)) {
  v <- vars[i]
  x <- min_neg[[v]]
  
  # numeric vs categorical branch
  if(is.numeric(x)) {
    res$p[i] <- t.test(x ~ min_neg$treat)$p.value
  } else {
    tbl <- table(x, min_neg$treat)
    res$p[i] <- chisq.test(tbl)$p.value
  }
}

print(res)

# Appendix B Regressions

# Texas Regressions

txdata <- read.csv("./tex.csv")
library(dplyr)
library(stargazer)
tx_pos <- txdata %>% filter(treatment < 2)
tx_neg <- txdata %>% filter(treatment > 1) %>% mutate(treatment = treatment - 2)

modelper_pos1 <- lm(data = tx_pos, personal_ben ~ treatment + Gender + Age + Mexican + Education +  Income + Born.in.U.S. + Partisanship + Texas.Id + American.Id + txknw + Travel + Family + SDO)
modelper_neg1 <- lm(data = tx_neg, personal_ben ~ treatment + Gender + Age + Mexican + Education +  Income + Born.in.U.S. + Partisanship + Texas.Id + American.Id + txknw + Travel + Family + SDO)
modeltx_pos2 <- lm(data = tx_pos, Tx_ben ~ treatment + Gender + Age + Mexican + Education +  Income + Born.in.U.S. + Partisanship + Texas.Id + American.Id + txknw + Travel + Family + SDO)
modeltx_neg2 <- lm(data = tx_neg, Tx_ben ~ treatment + Gender + Age + Mexican + Education +  Income + Born.in.U.S. + Partisanship + Texas.Id + American.Id + txknw + Travel + Family + SDO)
modelsup_pos3 <- lm(data = tx_pos, support ~ treatment + Gender + Age + Mexican + Education +  Income + Born.in.U.S. + Partisanship + Texas.Id + American.Id + txknw + Travel + Family + SDO)
modelsup_neg3 <- lm(data = tx_neg, support ~ treatment + Gender + Age + Mexican + Education +  Income + Born.in.U.S. + Partisanship + Texas.Id + American.Id + txknw + Travel + Family + SDO)
modelmed_pos4 <- lm(data = tx_pos, support ~ treatment + Gender + Age + Mexican + Education +  Income + Born.in.U.S. + Partisanship + Texas.Id + American.Id + txknw + Travel + Family + SDO + Tx_ben)
modelmed_neg4 <- lm(data = tx_neg, support ~ treatment + Gender + Age + Mexican + Education +  Income + Born.in.U.S. + Partisanship + Texas.Id + American.Id + txknw + Travel + Family + SDO + Tx_ben)
modelmed_pos5 <- lm(data = tx_pos, support ~ treatment + Gender + Age + Mexican + Education +  Income + Born.in.U.S. + Partisanship + Texas.Id + American.Id + txknw + Travel + Family + SDO + Tx_ben + U.S._ben)
modelmed_neg5 <- lm(data = tx_neg, support ~ treatment + Gender + Age + Mexican + Education +  Income + Born.in.U.S. + Partisanship + Texas.Id + American.Id + txknw + Travel + Family + SDO + Tx_ben + U.S._ben)

# For Study 1a Regression Table in Appendix B
stargazer(modelper_pos1,modeltx_pos2,modelsup_pos3, font.size = "footnotesize", title="Treatment Effects with Controls (Experiment 2)", dep.var.labels = c("Personal Benefit","Belief: Texas Benefit","Trade Support"), dep.var.caption = "", covariate.labels = c("Treatment", "Gender","Age","Mexican Ancestry (Self-Reported)","Education","Income","American Born","Partisanship","Texas Id Strength", "U.S. Id Strength", "Texas Knowledge/Socialization", "Frequency of Contact (Travel to Mexico)", "Family in Mexico", "SDO/16", "Belief: Texas Benefit (Mediator/Mechanism)"), single.row=F, omit=c("Constant"), omit.stat=c("ser","f"), no.space = F, column.sep.width = "5pt", style = "apsr", type = 'latex', header = FALSE, report=('vc*p'))

# For Study 2a Regression Table in Appendix B
stargazer(modelper_neg1,modeltx_neg2,modelsup_neg3, font.size = "footnotesize", title="Treatment Effects with Controls (Experiment 2)", dep.var.labels = c("Personal Benefit","Belief: Texas Benefit","Trade Support"), dep.var.caption = "", covariate.labels = c("Treatment", "Gender","Age","Mexican Ancestry (Self-Reported)","Education","Income","American Born","Partisanship","Texas Id Strength", "U.S. Id Strength", "Texas Knowledge/Socialization", "Frequency of Contact (Travel to Mexico)", "Family in Mexico", "SDO/16", "Belief: Texas Benefit (Mediator/Mechanism)"), single.row=F, omit=c("Constant"), omit.stat=c("ser","f"), no.space = F, column.sep.width = "5pt", style = "apsr", type = 'latex', header = FALSE, report=('vc*p'))

# Minnesota Regressions

mindata <- read.csv("./min.csv")
library(dplyr)
library(stargazer)
min_pos <- mindata %>% filter(treat < 3) %>% mutate(treat = treat - 1)
min_neg <- mindata %>% filter(treat > 2) %>% mutate(treat = treat - 3)

modelper_pos1 <- lm(data = min_pos, personal_ben ~ treat + Gender + Age + white + education +  income + us_born + partisanship + Minnesota.Id + American.Id + minkn + travel_abroad + family_abroad + fam_visit + SDO + twin_cities)
modelper_neg1 <- lm(data = min_neg, personal_ben ~ treat + Gender + Age + white + education +  income + us_born + partisanship + Minnesota.Id + American.Id + minkn + travel_abroad + family_abroad + fam_visit + SDO + twin_cities)
modelmin_pos2 <- lm(data = min_pos, min_ben ~ treat + Gender + Age + white + education +  income + us_born + partisanship + Minnesota.Id + American.Id + minkn + travel_abroad + family_abroad + fam_visit + SDO + twin_cities)
modelmin_neg2 <- lm(data = min_neg, min_ben ~ treat + Gender + Age + white + education +  income + us_born + partisanship + Minnesota.Id + American.Id + minkn + travel_abroad + family_abroad + fam_visit + SDO + twin_cities)
modelsup_pos3 <- lm(data = min_pos, support ~ treat + Gender + Age + white + education +  income + us_born + partisanship + Minnesota.Id + American.Id + minkn + travel_abroad + family_abroad + fam_visit + SDO + twin_cities)
modelsup_neg3 <- lm(data = min_neg, support ~ treat + Gender + Age + white + education +  income + us_born + partisanship + Minnesota.Id + American.Id + minkn + travel_abroad + family_abroad + fam_visit + SDO + twin_cities)

# For Study 1b Regression Table in Appendix B
stargazer(modelper_pos1,modelmin_pos2,modelsup_pos3, font.size = "footnotesize", title="Treatment Effects with Controls (Experiment 4)", dep.var.labels = c("Personal Benefit","Belief: Minn. Benefit","Trade Support"), dep.var.caption = "", covariate.labels = c("Treatment", "Gender","Age","White (Self-Reported)","Education","Income","American Born","Partisanship","Minnesota Id Strength", "U.S. Id Strength", "Minn. Knowledge/Socialization", "Frequency of Travel Abroad", "Family Abroad", "Frequency of Family Visits", "SDO/16", "Twin Cities Dummy (1/0)"), single.row=F, omit=c("Constant"), omit.stat=c("ser","f"), no.space = F, column.sep.width = "5pt", style = "apsr", type = 'latex', header = FALSE, report=('vc*p'))

# For Study 2b Regression Table in Appendix B
stargazer(modelper_neg1,modelmin_neg2,modelsup_neg3, font.size = "footnotesize", title="Treatment Effects with Controls (Experiment 4)", dep.var.labels = c("Personal Benefit","Belief: Minn. Benefit","Trade Support"), dep.var.caption = "", covariate.labels = c("Treatment", "Gender","Age","White (Self-Reported)","Education","Income","American Born","Partisanship","Minnesota Id Strength", "U.S. Id Strength", "Minn. Knowledge/Socialization", "Frequency of Travel Abroad", "Family Abroad", "Frequency of Family Visits", "SDO/16", "Twin Cities Dummy (1/0)"), single.row=F, omit=c("Constant"), omit.stat=c("ser","f"), no.space = F, column.sep.width = "5pt", style = "apsr", type = 'latex', header = FALSE, report=('vc*p'))

# Appendix C Summary Statistics

txdata <- read.csv("./tex.csv")
library(dplyr)
library(stargazer)
tx_pos <- txdata %>% filter(treatment < 2)
tx_neg <- txdata %>% filter(treatment > 1) %>% mutate(treatment = treatment - 2)

mindata <- read.csv("./min.csv")
min_neg <- mindata %>% filter(treat > 2) %>% mutate(treat = treat - 3)
min_pos <- mindata %>% filter(treat < 3) %>% mutate(treat = treat - 1)

# MISSING the original code for summary statistics, stargazer(dataframe name) will return the same results raw though.
# (variables will not be named neatly as in the original appendix, but this will not prevent one from matching them.)

stargazer(tx_pos)
stargazer(tx_neg)
stargazer(min_neg)
stargazer(min_pos)

# Appendix D

txdata <- read.csv("./tex.csv")
library(dplyr)
library(stargazer)
tx_pos <- txdata %>% filter(treatment < 2)
tx_neg <- txdata %>% filter(treatment > 1) %>% mutate(treatment = treatment - 2)

# numbers to recreate tables---tables and CIs input manually in LaTeX
summary(lm(data=tx_pos, support ~ treatment))
summary(lm(data=tx_pos, personal_ben ~ treatment))
summary(lm(data=tx_pos, Tx_ben ~ treatment))

summary(lm(data=tx_neg, support ~ treatment))
summary(lm(data=tx_neg, personal_ben ~ treatment))
summary(lm(data=tx_neg, Tx_ben ~ treatment))

mindata <- read.csv("./min.csv")
min_neg <- mindata %>% filter(treat > 2) %>% mutate(treat = treat - 3)
min_pos <- mindata %>% filter(treat < 3) %>% mutate(treat = treat - 1)

# numbers to recreate tables---tables and CIs input manually in LaTeX
summary(lm(min_ben ~ treat, data = min_pos))
summary(lm(support ~ treat, data = min_pos))
summary(lm(personal_ben ~ treat, data = min_pos))

summary(lm(min_ben ~ treat, data = min_neg))
summary(lm(support ~ treat, data = min_neg))
summary(lm(personal_ben ~ treat, data = min_neg))

# Appendix E
# The code below will recreate the results from both the main text Figure 5 and appendix E. B/c Appendix E
# includes the sensiticity analysis for the main text mediation analysis these are difficult to disentangle but all results are there.

txdata <- read.csv("./tex.csv")
library(dplyr)
library(stargazer)
tx_pos <- txdata %>% filter(treatment < 2)
tx_neg <- txdata %>% filter(treatment > 1) %>% mutate(treatment = treatment - 2)

library(mediation)
library(sandwich)
library(ggplot2)
mindata <- read.csv("./min.csv")
min_neg <- mindata %>% filter(treat > 2) %>% mutate(treat = treat - 3)
min_pos <- mindata %>% filter(treat < 3) %>% mutate(treat = treat - 1)
# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 5
set.seed(5)
med.fit.postex <- lm(Tx_ben ~ treatment + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw , data = tx_pos)
out.fit.postex <- lm(support ~ treatment + Tx_ben + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw, data = tx_pos)

med.out.postex <- mediate(med.fit.postex, out.fit.postex, treat = "treatment", boot = TRUE, mediator = "Tx_ben", robustSE = FALSE)
summary(med.out.postex)
sens.out1 <- medsens(med.out.postex, rho.by = 0.1, effect.type = "indirect")

for_graph1 <- data.frame(p = sens.out1$rho, avg_med_eff0 = sens.out1$d0, avg_med_eff0_lower = sens.out1$lower.d0, avg_med_eff0_upper = sens.out1$upper.d0, avg_med_eff1 = sens.out1$d1, avg_med_eff1_lower = sens.out1$lower.d1, avg_med_eff1_upper = sens.out1$upper.d1)

a <- ggplot(data = for_graph1, aes(x = p, y = avg_med_eff0)) + geom_line(col='blue') + geom_ribbon(aes(ymin = avg_med_eff0 - (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645), ymax =  avg_med_eff0 + (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645)), alpha = 0.1) + ylab("") + xlab("") + ggtitle("Study 1a") + geom_hline(yintercept=0, linetype="dashed", color = "red")

# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 2
set.seed(2)
med.fit.negtex <- lm(Tx_ben ~ treatment + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw, data = tx_neg)
out.fit.negtex <- lm(support ~ treatment + Tx_ben + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw, data = tx_neg)

med.out.negtex <- mediate(med.fit.negtex, out.fit.negtex, treat = "treatment", boot = TRUE, mediator = "Tx_ben", robustSE = FALSE)
summary(med.out.negtex)
sens.out2 <- medsens(med.out.negtex, rho.by = 0.1, effect.type = "indirect")

for_graph2 <- data.frame(p = sens.out2$rho, avg_med_eff0 = sens.out2$d0, avg_med_eff0_lower = sens.out2$lower.d0, avg_med_eff0_upper = sens.out2$upper.d0, avg_med_eff1 = sens.out2$d1, avg_med_eff1_lower = sens.out2$lower.d1, avg_med_eff1_upper = sens.out2$upper.d1)

b <- ggplot(data = for_graph2, aes(x = p, y = avg_med_eff0)) + geom_line(col='blue') + geom_ribbon(aes(ymin = avg_med_eff0 - (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645), ymax =  avg_med_eff0 + (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645)), alpha = 0.1) + ylab("") + xlab("") + ggtitle("Study 2a") + geom_hline(yintercept=0, linetype="dashed", color = "red")

# summary(med.out.pos)
# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 11
set.seed(11)
med.fit.posmin <- lm(min_ben ~ treat + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn, data = min_pos)
out.fit.posmin <- lm(support ~ treat + min_ben + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn, data = min_pos)

med.out.posmin <- mediate(med.fit.posmin, out.fit.posmin, treat = "treat", boot = TRUE, mediator = "min_ben", robustSE = FALSE)
summary(med.out.posmin)
sens.out3 <- medsens(med.out.posmin, rho.by = 0.1, effect.type = "indirect")

for_graph3 <- data.frame(p = sens.out3$rho, avg_med_eff0 = sens.out3$d0, avg_med_eff0_lower = sens.out3$lower.d0, avg_med_eff0_upper = sens.out3$upper.d0, avg_med_eff1 = sens.out3$d1, avg_med_eff1_lower = sens.out3$lower.d1, avg_med_eff1_upper = sens.out3$upper.d1)

c <- ggplot(data = for_graph3, aes(x = p, y = avg_med_eff0)) + geom_line(col='blue') + geom_ribbon(aes(ymin = avg_med_eff0 - (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645), ymax =  avg_med_eff0 + (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645)), alpha = 0.1) + ylab("") + xlab("") + ggtitle("Study 1b") + geom_hline(yintercept=0, linetype="dashed", color = "red")

# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 17
set.seed(17)
med.fit.negmin <- lm(min_ben ~ treat + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn, data = min_neg)
out.fit.negmin <- lm(support ~ treat + min_ben + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn, data = min_neg)

med.out.negmin <- mediate(med.fit.negmin, out.fit.negmin, treat = "treat", boot = TRUE, mediator = "min_ben", robustSE = FALSE)
summary(med.out.negmin)
sens.out4 <- medsens(med.out.negmin, rho.by = 0.1, effect.type = "indirect")

for_graph4 <- data.frame(p = sens.out4$rho, avg_med_eff0 = sens.out4$d0, avg_med_eff0_lower = sens.out4$lower.d0, avg_med_eff0_upper = sens.out4$upper.d0, avg_med_eff1 = sens.out4$d1, avg_med_eff1_lower = sens.out4$lower.d1, avg_med_eff1_upper = sens.out4$upper.d1)

d <- ggplot(data = for_graph4, aes(x = p, y = avg_med_eff0)) + geom_line(col='blue') + geom_ribbon(aes(ymin = avg_med_eff0 - (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645), ymax =  avg_med_eff0 + (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645)), alpha = 0.1) + ylab("") + xlab("") + ggtitle("Study 2b") + geom_hline(yintercept=0, linetype="dashed", color = "red")

par(mfrow = c(2, 4))
plot(sens.out1, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("negative"), levels = NULL, main = "Study 1a")

plot(sens.out3, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("negative"), levels = NULL, main = "Study 1b")

plot(sens.out2, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("negative"), levels = NULL, main = "Study 2a")

plot(sens.out4, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("negative"), levels = NULL, main = "Study 2b")

plot(sens.out1, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("positive"), levels = NULL, main = "Study 1a")

plot(sens.out3, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("positive"), levels = NULL, main = "Study 1b")

plot(sens.out2, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("positive"), levels = NULL, main = "Study 2a")

plot(sens.out4, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("positive"), levels = NULL, main = "Study 2b")

# Appendix F
# Appendix F is the exact same as Appendix E except controlling for perception of personal benefit.

txdata <- read.csv("./tex.csv")
library(dplyr)
library(stargazer)
tx_pos <- txdata %>% filter(treatment < 2)
tx_neg <- txdata %>% filter(treatment > 1) %>% mutate(treatment = treatment - 2)

library(mediation)
library(sandwich)
library(ggplot2)
mindata <- read.csv("./min.csv")
min_neg <- mindata %>% filter(treat > 2) %>% mutate(treat = treat - 3)
min_pos <- mindata %>% filter(treat < 3) %>% mutate(treat = treat - 1)
# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 5
set.seed(5)
med.fit.postex <- lm(Tx_ben ~ treatment + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw + personal_ben, data = tx_pos)
out.fit.postex <- lm(support ~ treatment + Tx_ben + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw + personal_ben, data = tx_pos)

med.out.postex <- mediate(med.fit.postex, out.fit.postex, treat = "treatment", boot = TRUE, mediator = "Tx_ben", robustSE = FALSE)
summary(med.out.postex)
sens.out1 <- medsens(med.out.postex, rho.by = 0.1, effect.type = "indirect")

for_graph1 <- data.frame(p = sens.out1$rho, avg_med_eff0 = sens.out1$d0, avg_med_eff0_lower = sens.out1$lower.d0, avg_med_eff0_upper = sens.out1$upper.d0, avg_med_eff1 = sens.out1$d1, avg_med_eff1_lower = sens.out1$lower.d1, avg_med_eff1_upper = sens.out1$upper.d1)

a <- ggplot(data = for_graph1, aes(x = p, y = avg_med_eff0)) + geom_line(col='blue') + geom_ribbon(aes(ymin = avg_med_eff0 - (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645), ymax =  avg_med_eff0 + (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645)), alpha = 0.1) + ylab("") + xlab("") + ggtitle("Study 1a") + geom_hline(yintercept=0, linetype="dashed", color = "red")

# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 2
set.seed(2)
med.fit.negtex <- lm(Tx_ben ~ treatment + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw + personal_ben, data = tx_neg)
out.fit.negtex <- lm(support ~ treatment + Tx_ben + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw + personal_ben, data = tx_neg)

med.out.negtex <- mediate(med.fit.negtex, out.fit.negtex, treat = "treatment", boot = TRUE, mediator = "Tx_ben", robustSE = FALSE)
summary(med.out.negtex)
sens.out2 <- medsens(med.out.negtex, rho.by = 0.1, effect.type = "indirect")

for_graph2 <- data.frame(p = sens.out2$rho, avg_med_eff0 = sens.out2$d0, avg_med_eff0_lower = sens.out2$lower.d0, avg_med_eff0_upper = sens.out2$upper.d0, avg_med_eff1 = sens.out2$d1, avg_med_eff1_lower = sens.out2$lower.d1, avg_med_eff1_upper = sens.out2$upper.d1)

b <- ggplot(data = for_graph2, aes(x = p, y = avg_med_eff0)) + geom_line(col='blue') + geom_ribbon(aes(ymin = avg_med_eff0 - (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645), ymax =  avg_med_eff0 + (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645)), alpha = 0.1) + ylab("") + xlab("") + ggtitle("Study 2a") + geom_hline(yintercept=0, linetype="dashed", color = "red")

# summary(med.out.pos)
# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 11
set.seed(11)
med.fit.posmin <- lm(min_ben ~ treat + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn + personal_ben, data = min_pos)
out.fit.posmin <- lm(support ~ treat + min_ben + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn + personal_ben, data = min_pos)

med.out.posmin <- mediate(med.fit.posmin, out.fit.posmin, treat = "treat", boot = TRUE, mediator = "min_ben", robustSE = FALSE)
summary(med.out.posmin)
sens.out3 <- medsens(med.out.posmin, rho.by = 0.1, effect.type = "indirect")

for_graph3 <- data.frame(p = sens.out3$rho, avg_med_eff0 = sens.out3$d0, avg_med_eff0_lower = sens.out3$lower.d0, avg_med_eff0_upper = sens.out3$upper.d0, avg_med_eff1 = sens.out3$d1, avg_med_eff1_lower = sens.out3$lower.d1, avg_med_eff1_upper = sens.out3$upper.d1)

c <- ggplot(data = for_graph3, aes(x = p, y = avg_med_eff0)) + geom_line(col='blue') + geom_ribbon(aes(ymin = avg_med_eff0 - (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645), ymax =  avg_med_eff0 + (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645)), alpha = 0.1) + ylab("") + xlab("") + ggtitle("Study 1b") + geom_hline(yintercept=0, linetype="dashed", color = "red")

# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 17
set.seed(17)
med.fit.negmin <- lm(min_ben ~ treat + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn + personal_ben, data = min_neg)
out.fit.negmin <- lm(support ~ treat + min_ben + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn + personal_ben, data = min_neg)

med.out.negmin <- mediate(med.fit.negmin, out.fit.negmin, treat = "treat", boot = TRUE, mediator = "min_ben", robustSE = FALSE)
summary(med.out.negmin)
sens.out4 <- medsens(med.out.negmin, rho.by = 0.1, effect.type = "indirect")

for_graph4 <- data.frame(p = sens.out4$rho, avg_med_eff0 = sens.out4$d0, avg_med_eff0_lower = sens.out4$lower.d0, avg_med_eff0_upper = sens.out4$upper.d0, avg_med_eff1 = sens.out4$d1, avg_med_eff1_lower = sens.out4$lower.d1, avg_med_eff1_upper = sens.out4$upper.d1)

d <- ggplot(data = for_graph4, aes(x = p, y = avg_med_eff0)) + geom_line(col='blue') + geom_ribbon(aes(ymin = avg_med_eff0 - (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645), ymax =  avg_med_eff0 + (abs(avg_med_eff0 - avg_med_eff0_lower)/1.96*1.645)), alpha = 0.1) + ylab("") + xlab("") + ggtitle("Study 2b") + geom_hline(yintercept=0, linetype="dashed", color = "red")

par(mfrow = c(2, 4))
plot(sens.out1, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("negative"), levels = NULL, main = "Study 1a")

plot(sens.out3, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("negative"), levels = NULL, main = "Study 1b")

plot(sens.out2, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("negative"), levels = NULL, main = "Study 2a")

plot(sens.out4, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("negative"), levels = NULL, main = "Study 2b")

plot(sens.out1, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("positive"), levels = NULL, main = "Study 1a")

plot(sens.out3, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("positive"), levels = NULL, main = "Study 1b")

plot(sens.out2, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("positive"), levels = NULL, main = "Study 2a")

plot(sens.out4, sens.par = c("R2"),
     r.type = c("residual"), sign.prod = c("positive"), levels = NULL, main = "Study 2b")

# Appendix G

txdata <- read.csv("./tex.csv")
library(dplyr)
library(stargazer)
tx_pos <- txdata %>% filter(treatment < 2)
tx_neg <- txdata %>% filter(treatment > 1) %>% mutate(treatment = treatment - 2)

# numbers to recreate figures--pictures made manually in excel
summary(lm(data=tx_pos, support ~ treatment))
summary(lm(data=tx_pos, Tx_ben ~ treatment))
summary(lm(data=tx_pos, support + Tx_ben ~ treatment))

summary(lm(data=tx_neg, support ~ treatment))
summary(lm(data=tx_neg, Tx_ben ~ treatment))
summary(lm(data=tx_neg, support + Tx_ben ~ treatment))

mindata <- read.csv("./min.csv")
min_neg <- mindata %>% filter(treat > 2) %>% mutate(treat = treat - 3)
min_pos <- mindata %>% filter(treat < 3) %>% mutate(treat = treat - 1)

# numbers to recreate figures--pictures made manually in excel
summary(lm(min_ben ~ treat, data = min_pos))
summary(lm(support ~ treat, data = min_pos))
summary(lm(min_ben + support ~ treat, data = min_pos))

summary(lm(min_ben ~ treat, data = min_neg))
summary(lm(support ~ treat, data = min_neg))
summary(lm(min_ben + support ~ treat, data = min_neg))

# Appendix I

txdata <- read.csv("./tex.csv")
library(dplyr)
library(stargazer)
tx_pos <- txdata %>% filter(treatment < 2)
tx_neg <- txdata %>% filter(treatment > 1) %>% mutate(treatment = treatment - 2)

txdatatable2 <- lm(data = txdata, support ~ Gender + Age + Mexican + Education +  Income + Born.in.U.S. + Partisanship + American.Id + SDO)
mindatatable3 <- lm(data = mindata, support ~ Gender + Age + white + education + income + us_born + partisanship + American.Id + SDO)

stargazer(txdatatable2,mindatatable3)

